Semiclassical approach to Bose-Einstein condensates in a triple well potential 
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We present a new approach for the analysis of Bose-Einstein condensates in a few mode approxi- 
mation. This method has already been used to successfully analyze the vibrational modes in various 
molecular systems and offers a new perspective on the dynamics in many particle bosonic systems. 
We discuss a system consisting of a Bose-Einstein condensate in a triple well potential. Such sys- 
tems correspond to classical Hamiltonian systems with three degrees of freedom. The semiclassical 
approach allows a simple visualization of the eigenstates of the quantum system referring to the 
underlying classical dynamics. From this classification we can read off the dynamical properties of 
the eigenstates such as particle exchange between the wells and entanglement without further calcu- 
lations. In addition, this approach offers new insights into the validity of the mean-field description 
of the many particle system by the Gross-Pitaevskii equation, since we make use of exactly this 
correspondence in our semiclassical analysis. We choose a three mode system in order to visualize 
it easily and, moreover, to have a sufficiently interesting structure, although the method can also 
be extended to higher dimensional systems. 

PACS numbers: 03.75.Kk, 03.75.Lm, 03.65.Sq 



I. INTRODUCTION 

Bose-Einstein condensates form one of the main topics 
of research at the moment. One reason for this enor- 
mous interest is the fact that they combine concepts and 
techniques from different areas of physics, such as quan- 
tum optics, condensed matter physics, molecular physics 
and quantum chaos. On the experimental side, there has 
been a remarkable progress in confining and manipulat- 
ing Bose-Einstein condensates [1, 2], which has stimu- 
lated the theoretical research in the area. 

There has been a large number of previous studies an- 
alyzing the dynamics of Bose-Einstein condensates in a 
double well potential using a mean-field approach, the 
Gross-Pitaevskii equation. In [3] Smerzi et al. discussed 
the occurrence of macroscopic self-trapping within one 
well. The behavior of the system during an adiabatic 
change of parameters was studied in [4-6] and general- 
izations of the linear two level crossing scenarios and the 
Landau-Zener formula were analyzed. Another line of 
investigation considers the mesoscopic regime in which 
the quantum and the classical, i.e. mean-field descrip- 
tions overlap and therefore semiclassical techniques can 
be used to study the system [7-12]. 

In this paper we present a semiclassical technique to 
analyze the spectral properties of a Bose-Einstein con- 
densate. This method has already been used to describe 
the vibrational spectra of molecules [13-16] and provides 
an intuitive picture for the at first sight uninterpretable 
spectra. In this sense, the dynamics of a Bose-Einstein 
condensate in a triple well potential is analogous to the 
vibrations of a tri-atomic molecule. Using the geometri- 
cal language of classical mechanics to describe the quan- 



tum system, we introduce a simple semiclassical method 
of visualization and classification of the quantum eigen- 
states which allows a characterization of the dynamics 
of the system. Furthermore, we investigate how far the 
correspondence between the mean-field system and the 
quantum many body system can be extended when the 
number of particles decreases. 

For our studies we consider a Bose-Einstein condensate 
in a triple well potential, since the technique easily allows 
us to go beyond the standard double well potential anal- 
ysis. A triple well potential has a much richer structure 
[11, 17-21] and the power of the method can be shown 
without loss of clarity, still allowing a direct visualization 
of all relevant structures. 



II. THE MODEL 

In the following analysis we consider a system consist- 
ing of bosonic particles in an external periodic potential 
V{r) = V{r^rf) withrY= hdiei+l2d2e2+hd3e3, h e¥i 
and dk G R. If a weak two-particle point-like interaction 
is assumed, then the Hamiltonian in second quantization 
can be written as 



H 



d^r^^r) 



2m 



A + V{r) $(f) 
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Here, m is the particle mass, g ~ Anasfi^ /m is the cou- 
pling constant describing two-body interactions and as is 
the s-wave scattering length. For a repulsive interaction, 
g is positive while for an attractive interaction g takes a 
negative value. For the rest of the paper we choose scaled 
units with h — m — 1. The field operator <&(r) can be 



expanded in terms of bosonic annihilation operators, 

where we assume that the basis functions {(j)n,m} of the 
one-particle Hilbert space are exponentially localized in 
space and real, as is the case for the Wannier functions 
[22]. The index n describes basis functions in different 
wells and we will take into account only three different 
wells in order to model the three well potential. The 
second index m labels the excited states within a sin- 
gle well. Assuming Bose-Einstcin condensates, we can 
restrict ourselves to the lowest energy state m — 1 and 
neglect higher excited states (see also [7] for a careful 
discussion of this topic for a two well potential). Ex- 
perimentally such a system was realized in [1] for a two 
well potential but the technique can in principle also be 
extended to three wells. 

Expanding the Hamiltonian in this basis and neglect- 
ing fourth order terms in the creation and annihila- 
tion operators from different basis functions (modes) 
yields the well-known Bose-Hubbard Hamiltonian [23] re- 
stricted to three wells. So, the Hamiltonian can be writ- 
ten in a symmetrized form as 



H = Ho + W 



(3) 



with 



- a|ai + aiOi 03^2 + 0,202 

n-a = i^l h ^2 

+ ^3 ^ +x,^ ^ 

^^ /4a2 + a2ah%^j4«3 + a34y ^4^ 

W = — {a{a2 + ai,ai) — [a'^a'i + a!^a2) ■ (5) 

Here we neglect a constant energy shift. For conve- 
nience, we will choose the nonlinear interaction strengths 
Xj equal for each well in the following sections which is 
also in accordance with experimental realizations. Such 
Hamiltonians have already been studied in great de- 
tail for the more restrictive two mode model (e.g. in 
[10, 12, 24]). 

The Hamiltonian commutes with the particle number 
operator N ~ ni+n2 +^.3 which expresses the conserva- 
tion of the total number of particles. The symmetrized 
form is more convenient when considering the semiclas- 
sical limit, as will become clear in the next paragraph. 
Hamiltonians of this kind have been used in molecular 
physics in order to describe and assign vibrational spec- 
tra [13, 16]. In the molecular case they describe all kinds 
of vibrational degrees of freedom like stretches, bends, 
torsions etc. and include various resonant interactions 
corresponding to different simple rational ratios between 
the frequencies. The conserved particle number in our 



case of Eq. (3) corresponds to the polyad-type conserved 
quantities in the molecular systems. 

For the case of 30 particles considered in the following, 
it is an easy numerical task to diagonalize the Hamilto- 
nian matrix and thus solve the problem. However, one 
cannot understand the underlying structure of this sys- 
tem from numerical values alone. The aim of this paper 
is to present a method which allows an easy visual char- 
acterization of the eigenvectors of the Hamiltonian, using 
the close correspondence with the classical system. 



A. The classical system 

Essential for our semiclassical classification and assign- 
ment of quantum states is a comparison between the 
quantum states and the corresponding classical dynam- 
ics. To this end, the first step is the construction of the 
classical Hamiltonian function, which corresponds to the 
quantum Hamiltonian given in Eqs. (3)-(5). This is done 
by Heisenberg's substitution rules [25] 



flfc 



he'' 



he' 



(6) 



There are two different lines of argumentation for this 
substitution. First, it is exact for the harmonic oscil- 
lator where the well known classical Hamiltonian ul is 
obtained by the replacement of the symmetrized product 
of an annihilation and a creation operator by the classical 
action. This implies the correspondence 



(7) 



between the classical action / and the quantum num- 
ber n of the oscillator (/ is here measured in units of 
h). This correspondence of Eq. (7) is also a result of the 
application of the semiclassical Bohr-Sommerfeld quan- 
tization rules to the harmonic oscillator. In more general 
cases we have to generalize the Bohr-Sommcrfeld method 
to the EBK quantization. Then the argument holds for 
any bound system of any number of degrees of freedom 
as long as the system is close to integrable (for general 
background information on semiclassics see [26]). In gen- 
eral, the semiclassical methods give results correct in the 
lowest two orders in h (orders and 1) and cause errors 
of order h"^. The application of the substitution rules of 
Eq. (6) to the quantum Hamiltonian of Eqs. (3)-(5) gives 

H(Lpi,(f2,f3,Il,l2,h) 

= Hoih, I2, h) + W(lpi, (f2, ip3, h, h, h) 

= ujili+ L02 I2 +UJ3I3+ xi if + X2II + X3 ll (8) 

- ki2\/lll2 COs((pi - IP2) - k2?.\/hh COs((y92 - 1^3.) ■ 

A Hamiltonian for the same system but expanded in an- 
other basis was analyzed in [11]. This function can be in- 
terpreted as the Hamiltonian of a classical system of three 
coupled anharmonic oscillators described in action-angle 



variables (ft G [0, 27r) and Ik > 0, where fc = 1, 2, 3. As a 
niethod to construct the corresponding classical Haniil- 
tonian, the substitution rules of Eq. (6) always give the 
correct result since in this direction (quantum — > clas- 
sical) the correspondence is unique whenever it exist at 
all, in contrast to the other direction (classical -^ quan- 
tum) with its notorious h^ problems. At high excitation 
(large quantum numbers) there is a second argument for 
the semiclassical correspondence. The application of a 
creation or annihilation operator to a number state \n) 
has the effect 



a\n) ^ ^/n\n — 1)^ a) \n) = \/n + 1 |n + 1) 



(9) 



In the limit of a large quantum number n, the difference 
between n and n -I- 1 or n — 1 is irrelevant in the square 
roots as well as in the states and the operators can sim- 
ply be replaced by multiplication with the number ^/n. 
This argument holds for condensates where a large num- 
ber of particles goes into a superfluid state which is well 
described by a mean- field limit. This is in line with the 
standard argument of semiclassical behavior in the limit 
of large quantum numbers. Interestingly, for systems of 
coupled anharmonic oscillators the semiclassical treat- 
ment is very good also for low excitation numbers. In this 
limit, we approach the integrable harmonic limit where 
the Bohr-Sommerfeld treatment gives the correct result. 
The experience with molecular systems of the structure 
of Eqs. (3)-(5) shows that a semiclassical treatment of 
such systems is globally quite good in most cases. 

Accordingly, we base our method of semiclassical as- 
signment on this argument. Semiclassical arguments will 
be used later first to convert the eigenstates of the many- 
body Hamiltonian (3) into wave functions on the toroidal 
configuration space and second to compare these func- 
tions with important structures seen in the classical dy- 
namics. 

The integrable part _ffo of the Hamiltonian, which does 
not contain interactions between the three oscillators, 
leaves all actions unchanged. In contrast, W changes 
the values of the actions (particles in the wells) because 
of its dependence on angles and introduces interactions 
between the three oscillators. In this sense we call in the 
following W the interaction part of the Hamiltonian. In 
the picture of particles in the triple well, W describes 
tunneling terms between the various wells. 

The Poisson bracket between H and the observable 



K = h+l2+h 



(10) 



the total action, is equal to zero, which corresponds to 
the quantum mechanically conserved number of particles. 
Note that the numerical value of K differs by 3 • 1/2 from 
the value of N because of the zero point actions. The 
symmetry {_ff, K} ~ can be used to reduce the number 
of degrees of freedom from three to two by a canonical 
transformation. Using the generating function 

= Jl{<Pl - (p2) + J2{V3 - V2) + KlP2 (11) 



of the old angles {ipi,<f2,f'i) and the new actions 
( Ji, J2, K) results in the transformations (together with 
Eq. (10)) 

/l = Jl , /3 = J2 , (12) 

where (^1,1/12,^) are the new angles conjugate to 
iJi,J2,K).^ 

The Hamiltonian in the new coordinates is given by 

H = LJi Ji + LO2 [K - Ji - J2) + bJz J2 



Xi Jl +X2{K- Ji- J2f + 2:3 J^ 



ki2\/ Ji{K - Ji - J2) costAi 



- fc23 \J -hiK - Jl - J2) cos -4)2 , 
with corresponding equations of motions 
^1 = {uji + 2xiJi) - {uj2 + 2x2{K - Jl- J2)) 



(13) 
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fc23 
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K-J1-J2 



Jl 



Jl 



K-J1-J2 



■h 



■ cos 1/12 



K-JI-J2 

1P2 = (w3 + 20:3^2) - i0J2 + 2X2{K - Jl - J2)) 
_ ^23 

2 



K-JI-J2 



J2 



J2 



K-JI-J2 



fcl2 

2 



Jl 



K-J1-J2 



COSl/ii 



Jl = -^12 \/ Ji{K - Jl - J2) sinV'i , 



J2 = -^23 \/J2{K- Jl - J2) sinV'2 • 



cosi/'i 
(14) 

cos 1/^2 

(15) 

(16) 
(17) 



The classical configuration space is a two dimensional 
torus spanned by the two angles ipi and '02 • In order to 
compare the classical and the quantum system we have 
to represent the states as wave functions on the classical 
configuration space. The way to do this will be described 
in the following section. 



B. The quantum mechanical configuration space 

The angle variables can be introduced in the quantum 
system by using the set of functions 



1^1,^2,^3)= Ve*("^^^+"^^^+"^^^Mm,n2,n3), (18) 



E 

ni,n2,n3>0 



first introduced in molecular spectroscopy by Sibert and 
McCoy [13]. These functions are similar to the Bargmann 
states studied in [10] in the context of a Bose-Einstein 
condensate. This relation is well-known from the context 
of infinite lattices. There, the sum is taken from —00 to 
00 and corresponds to the representation of Bloch func- 
tions in terms of Wannier functions. The angle variables 



(/3i, (/?2 and (p3 span the Brillouin zone. However, in this 
example these functions arc not orthogonal due to the 
fact that for fixed N the sum is finite, 

— y^ g-«(«i(¥'i-'p'i)+"2(¥'2-y2)+"3(v3-v3)) ^ ng\ 

ni+n2+n3=N 

For a large particle number N the scalar product con- 
verges to a delta-comb. There is a considerable deviation 
for the value oi N = 30, which can play an important role 
when matrix elements are calculated. But here we use 
these functions only for visualization and not for further 
algebraic manipulations. The eigenfunctions of (3) have 
the form 



1$) 



ni+n2+n3=N 



ii,n2,n3 1^1,112,713) . 



(20) 



The coefficients Cni.n2,n3 can be obtained by a numeri- 
cal diagonalization in the number basis |rii,n2,n3). The 
eigenstates in the angle representation (18), i.e. the wave 
functions, are given by a Fourier series: 



{fl,f2,f3\'^) = Yl 



Jiniipi+n2(p2+n3(p3) 



ni+n2+n3=N 



(21) 

Finally, we can reduce the number of degrees of freedom 
in this representation by using the same coordinate trans- 
formation as in the classical case of Eq. (12). This leads 
to the expression 



$(V'l,V'2) = (V'l,^2|$) 



„iJVi9 



/ ^ (^7ii.N — ni—7i3,n3 ^ 
ni-\-n3<N 



i(niipi+n3'ip2) 



(22) 



The global phase factor e*^"* can be ignored in the follow- 
ing considerations. It must be emphasized that the sum 
includes only a finite number of terms due to the finite 
number of combinations of numbers rii, n2 and 71,3 which 
sum to N. Therefore the Fourier expansion in Eq. (22) 
has only a finite resolution. For a very small value of N, 
this sum has just a few terms, so that only the coarse 
grain structure can be explored; accordingly, the eigen- 
functions ^(ipi,ip2) show only diffuse structures. In our 
example, the configuration space of the reduced system 
is the two dimensional torus T^ with total volume 47r^. 
The total number of basis states for a given number of 
particles N is L = {N + 1){N + 2)/2. Accordingly, the 
eigenfunctions which are linear combinations of the L ba- 
sis functions can only show patterns with a resolution of 
the order Air"^ /L in the area or a resolution of the order 
2tt/N in each direction. With iV = 30 we have L = 496 
eigenstates, giving a resolution of approximately O.OTtt in 
each direction. 

The reinterpretation of the expansion of an eigen- 
state into number states as a Fourier series on the 



toroidal configuration space has the following semiclas- 
sical interpretation, where we write for the moment 
h explicitly into the equations: If one naively quan- 
tizes the classical canonical variables {(pk,Ik) using the 
Schrodinger quantization [(pk,Ii] — itiSki, which imposes 
Ik = {ih)^^d/{d(f>k), then functions f{(f>i,(f>2,^3) can 
be interpreted as wave functions in coordinate space. 
Of course, the Schrodinger quantization is correct only 
in Cartesian coordinates and it does not commute with 
canonical transformations, in general yielding errors of 
the order of fi?. Therefore the results have to be inter- 
preted semiclassically Note that due to our symmetric 
introduction of the quantum-classical correspondence in 
Eq. (7), the errors to first order in h cancel identically. 
Because of these considerations we call the wave function 
from Eq. (22) the semiclassical wave function. 

In many semiclassical investigations, Husimi functions 
are used to relate quantum wave functions of eigenstates 
to structures in the classical phase space. This is the ap- 
propriate and natural procedure if the usual position and 
momentum coordinates are used. It is less clear and in 
addition not necessary in our case where the whole dy- 
namics is treated in action-angle variables. Let us explain 
this point in some detail: The description of the system 
by a Hamiltonian of the functional structure of Eqs. (3)- 
(5) in the quantum case or Eq. (8) in the classical case 
only makes sense for bound systems, it is not appropriate 
to describe scattering systems. Therefore we restrict the 
following discussion to bound states only. For any bound 
eigenstate in the standard position space, there must be 
the same amount of wave running in one direction and in 
the opposite direction, otherwise it would not be a bound 
stationary state. Accordingly, the wave function can be 
chosen real. The phases of the wave function do not play 
any important role and do not help for the classifica- 
tion of the states. The canonically conjugate momenta 
have continuous values and Wigner or Husimi functions 
are defined without any problem on the classical phase 
space and indicate in many cases to which structure in 
the classical phase space some particular quantum state 
belongs. 

The situation is very different in action angle variables. 
Here the configuration space is a torus with its very dif- 
ferent global topology. This causes great difficulties to 
define the usual Wigner or Husimi functions. Because 
of the periodicity of the configuration coordinates, the 
corresponding canonically conjugate variables (here the 
actions) only have discrete values in the quantum dynam- 
ics. This makes it very tricky to convert the wave func- 
tion into something defined on the continuous classical 
phase space. On the other hand, we do not really need 
to do this, since we have the following simpler method 
to squeeze out of the wave functions information on the 
classical actions. Waves propagating in one direction on 
a torus always return to the starting point. Accordingly, 
wave functions for a bound state can have - and in fact 
in most cases do have - strong running wave contribu- 
tions and the phase of the function is essential and will 



be analyzed to help in the classification of the state. In 
a semiclassical spirit the phase of a wave function can be 
interpreted as a classical action integral and accordingly 
the gradient of the phase function gives the value of the 
canonically conjugate momentum which in this case is 
the action. If there is a sufficiently large patch of con- 
figuration space where the phase function comes close 
to a plane wave, then its gradient indicates the value of 
the actions which is represented by this part of the wave 
function. This provides a kind of lift of the wave function 
from configuration space into phase space. If there are 
closed loops on the torus along which the phase function 
is very regular (and this usually happens along density 
crests which run along the classical organizing center as 
will be explained in detail in section IV) then we interpret 
this as representing a motion of almost constant action 
along this loop. This idea is used to get longitudinal 
quantum numbers introduced in section IV. 



III. CLASSICAL DYNAMICS AND COUPLING 

SCHEMES 

Before we relate individual quantum states to guiding 
centers of the classical dynamics, we must get an overview 
of the classical dynamics and its skeleton. As an exam- 
ple, we discuss the classical dynamics for N = 30, i.e. for 
the value 31.5 of the classically conserved total action K. 
In the following, we choose parameter values lui = —lut, = 
0.1, UJ2 = 0, xi = a;2 = xa = 0.1 and fci^2 = fe,3 = 0.5, 
which lead to a quantum mechanical energy interval of 
[23.907,96.393]. The classical reduced system exists in 
the energy interval [22.476,99.1]. Furthermore, we mea- 
sure all energies with respect to the quantum mechanical 
ground state of iJ in Eq. (3), i.e. we subtract the quantum 
mechanical zero point Ho(l/2, 1/2, 1/2) — 0.075 from the 
classical energies in order to facilitate the comparison 
between classical and quantum dynamics. To represent 
the classical dynamics graphically, we show Poincare sec- 
tions in planes tpi — ^ with positive orientation i/ii > 0. 
If an initial condition (i/;2, J2) is chosen in the Poincare 
section, then we first have to reconstruct the four corre- 
sponding coordinates in the phase space in order to start 
a trajectory of the flow through this point. The two co- 
ordinates "02 and J2 coincide with the given coordinates 
in the domain of the Poincare map. The coordinate ipi is 
obtained by the intersection condition and the remaining 
coordinate Ji is calculated by an inversion of the Hamil- 
tonian function (13) with respect to the coordinate Ji for 
a fixed value of the energy and for the known values of 
the other three coordinates. Here some care is necessary 
since this inverse function is multivalued. First we fix one 
orientation of the domain, i.e. we always search for solu- 
tions with (hpi/dt > 0. In principle there can be several 
solutions with the same orientation and then it is neces- 
sary to ensure that all initial points used belong to the 
same branch. Poincare sections in planes 1/12 — constant 
look very similar to the ones in planes ■01 = constant. 



Therefore it is sufficient to restrict ourselves to sections 
in 01 = only. 

If the whole dynamics were governed by _ffo, then all 
actions would be constants of motion and all Poincare 
sections would be foliated by invariant lines J2 = 
constant. Including the interaction W between the wells 
(modes) into the dynamics has the following effects. In 
regions of the phase space, where none of the resonances 
contained in W has an important effect, the dynamics 
is in the KAM regime (see the extensive discussion of 
soft chaos in chapter 9 of [27]) and a large fraction of 
the phase space volume is still filled by invariant lines, 
which are continuous deformations of the invariant sur- 
faces J = constant of the unperturbed Hq dynamics. We 
call such invariant surfaces primary tori. This happens 
mainly in regions of phase space where the effective fre- 
quencies 



off _ ^-^0 



(23) 



are far from simple rational ratios, for which there is a 
corresponding resonance coupling in W ^ as explained in 
the next paragraph. For our particular choice of coupling 
terms in W, only 1:1 resonances are relevant. 

The effect of the coupling terms between the different 
modes can be described in the following way. Each term 
contains a cosine function whose argument is a difference 
between angles of the original degrees of freedom or one 
angle of the reduced system, see Eqs. (8) and (13). Be- 
cause in our special case the arguments are differences 
of two angles with the same weight, we say that these 
terms describe 1:1 resonant interactions between the two 
degrees of freedom. The right hand sides of the Hamil- 
tonian equations of motion (14) and (15) for the angles 
ij^k (fc = 1, 2) of the reduced system. 



dt 



dJk 



dW 
dJk 



(24) 



contain two contributions. The first consists of the dif- 
ference of two effective frequencies from Eq. (23), and 
the second is the derivative of the coupling terms with 
respect to the action, which contains cosine functions. 
First, let us assume that we change some parameter, e.g. 
fci,27 to see how coupling sets in. Further we assume that 
the difference between the effective frequencies, i.e. the 
angle independent term on the right hand side, is differ- 
ent from zero. Let us say it has the value v ^ Q. For 
a small value of fci.2 the angle dependent terms are not 
able to cancel v regardless of the value of the angles. The 
angle dependent terms have the maximal absolute value 
for angle values and ■n because of the dependence on co- 
sine functions. When fci_2 increases, then at one point it 
reaches a value, where the angle dependent terms are just 
able to cancel v. Then the angle ■0fc of the reduced sys- 
tem stops, ipk{t) = constant, and we call this frequency 
locking. This necessarily happens for angle values where 
the cosine functions have maximal absolute value, i.e. 



where the angles are or n. Whether the appropriate 
angle values are or tt depends on the signs of v and 
of the terms in front of the cosine functions. When the 
value of fci.2 is further increased, then there is a whole 
interval of angle values where locking is possible. The ac- 
tual dynamics of the locked motion then performs small 
oscillations around the angle values or tt. This will be 
seen in the numerical results of the classical dynamics. In 
the quantum dynamics the fluctuations around the cou- 
pling point of the angles are quantized and give rise to a 
discrete set of transversal quantum numbers, see section 
IV. 

If only one of these resonant couplings is strong, then 
the dynamics is still close to integrable, and a large part 
of the phase space volume is filled by invariant tori, which 
show up as invariant lines in the Poincare sections. How- 
ever, due to the rearrangement of phase space structures 
by the resonant coupling, the invariant surfaces in phase 
space are no longer primary tori, i.e. are no longer con- 
tinuous deformations of invariant surfaces of the Hq dy- 
namics. Large bundles of secondary tori appear which 
are organized around periodic orbits (in this case sta- 
ble, elliptic) representing the guiding centers for the new 
nonlinear modes. There are also corresponding unstable 
periodic orbits, which in the integrable case are repre- 
sented by separatrix crossings in Poincare sections. In 
the nonintegrable cases, the separatrices break and turn 
into homoclinic tangles, which become the central struc- 
tures of chaotic strips. However, if only one resonant 
coupling has a strong effect and the others are not im- 
portant, then the chaos strips are very thin and they still 
appear almost like separatrices. 

If two or more linearly independent resonant couplings 
arc strong, then chaos on large scales can appear. These 
regions in phase space are resonance overlap zones [28]. 
However, also in strongly chaotic regions of phase space 
there are still simple short periodic orbits (in this case 
unstable, normal hyperbolic or inverse hyperbolic) which 
act as guiding centers of the flow. Then the dynamics 
is chaotic but nevertheless the flow follows some guid- 
ing center on the average. This average flow is relevant 
for the comparison with quantum dynamics. Thus, also 
in the classically chaotic case we may find surprisingly 
simple and clean structures in a large part of the quan- 
tum wave functions. In such cases it can be appropriate 
to imagine simple idealized classical guiding centers and 
interpret the quantum states as quantum excitations of 
these idealized structures. 

Let us give a short estimate of the size of structures 
which are relevant for our semiclassical considerations. 
The range of action values is limited between and K 
due to Eq. (10), the angle can vary over an interval of 
length 27r. For each particular plot only a part of this 
range is energetically accessible in reality. Accordingly 
the size of the Poincare section is limited by 2ttK. For 
semiclassical investigations structures of a size of h or 
larger are relevant. We always use units in which h has 
the numerical value 1 and also the values of all actions 



should be interpreted as being given in units of h. There- 
fore structures in our Poincare plots are of interest in the 
following, if their size is at least in the order of one unit 
of action or has a relative size of l/K compared to the 
size of the maximally possible domain of the map. 

We perform almost all our calculations for the reduced 
system. On the other hand, the real object of interest is 
the original system of particles in three wells. Therefore 
we need a fast and easy method to transfer statements 
about the reduced system into the corresponding state- 
ments about the original system. We have called this 
procedure the lift in the previous work on molecular sys- 
tems [14-16]. Let us assume a trajectory of the reduced 
system is given and we want to reconstruct the corre- 
sponding trajectory of the original system. The first step 
of the procedure is the reconstruction of the cyclic angle. 
It is done rigorously by using the Hamiltonian equation 
of motion 

^ - ^ (25) 

dt~ dK- ^^^> 

The right hand side of this equation does not depend on 
?? but only on the known values of the other coordinates 
as function of time. Accordingly we get 'd{t) by a simple 
integration with respect to time. The experience with 
the molecular systems has shown that normally it is suf- 
ficient to approximate 'd{t) by t times a constant effective 
frequency. In our case i? is the only fast variable of the 
whole system and describes a fast oscillation superim- 
posed on the motion of the whole system. The initial 
value ^(0) is rather irrelevant. In contrast, the variables 
of the reduced system are slow variables describing the 
relative motion between the various degrees of freedom of 
the original system. The next step of the lift procedure 
is to undo the canonical transformation and to go back 
to the coordinates of the original system. In this second 
step the advantage of choosing the new actions equal to 
some of the old actions becomes evident. The knowledge 
of the actions in the reduced system and of the constant 
value of K gives immediately the values of the old actions, 
i.e. the values of the particle numbers in the three wells. 
Because of this simple connection between the actions of 
the reduced system and the actions of the original system 
we will switch very freely between the reduced and the 
original system in the following considerations. 

In our case we have in the interaction part W of the 
Hamiltonian 1:1 couplings between the degrees of free- 
dom 1 and 2 and between 2 and 3 respectively. Indi- 
rectly this also implies a 1:1 coupHng between the degrees 
of freedom 1 and 3. Accordingly, we have the following 
coupling schemes: 

Type (A) : If the effective frequencies are not very close 
to each other, then no interaction term can cause fre- 
quency and phase coupling, and all three modes run in- 
dependently with their own effective frequency. This is 
the KAM regime with many primary tori, where the mo- 
tion is of quasiperiodic type with three independent fre- 
quencies. The organization center of the reduced sys- 
tem is the complete configuration space T^. In Poincare 




FIG. 1: The classical reduced system for energy E — 55. 
(a) Poincare plot in the plane ■i/'i = for variables ip2 and 
J2 with Ji fixed by energy conservation, (b) Trajectory in 
a primary torus in the lower region of (a) for initial values 
{ipi,ip2, J2) = (0, n, 7.5). (c) Trajectory in a primary torus in 
the upper region of (a) for initial values (0, 0, 20.6). The points 
of the trajectories are given in equidistant time intervals At = 
0.01 in order to indicate the velocity by the distance between 
neighboring points. 



plots, we see many invariant lines which are continuous 
deformations of horizontal lines J2 — constant, i.e. of 
the invariant lines belonging to Hq . This type of motion 
appears mainly in the middle of the accessible energy 
interval for a given particle number. In Fig. 1 we give 
some numerical results for the energy E = 55. Part (a) 
shows the Poincare section and parts (b) and (c) show 
two segments of trajectories in the reduced configuration 
space. The domain of the Poincare map in (a) consists of 
two parts. The range of J2 values between approximately 
10.2 and 20.5 is not accessible at this energy. At values of 
J2 around 9, we see many primary tori. A segment (five 
revolutions in direction of ipi) of a typical trajectory be- 
longing to one of them is shown in part (b) of the figure. 
In the long run, the trajectory fills the whole configura- 
tion space quasiperiodically. In these primary tori the 
action J2 is smaller than the action Ji, so that the tra- 
jectories move faster in -01 direction than in 1/12 direction. 
The opposite happens on the primary tori lying around 
J2 values of 21. Here the J2 action is largest and there- 
fore the quasiperiodic trajectories run with higher speed 
in the V'2 direction. (For a numerical example, see a tra- 
jectory segment in Fig. 1(c)). The other structures seen 
in Fig. 1(a) belong to other types of motion, discussed 
below. 



Type (B) : If the effective frequencies of modes 2 and 
3 are close but that of the first mode is not close, then 
we expect that modes 2 and 3 are locked but mode 1 
is independent. The motion is then quasiperiodic with 
two independent frequencies. The organization center 
in the reduced system is a one dimensional curve with 
■02 = constant. In Poincare plots in the plane "01 = 0, we 
see secondary islands. This motion appears mainly for 
high energies. Figure 2 gives some numerical results for 
E = 80. Part (a) shows a Poincare section, again in the 
plane 0i ~ 0, and part (b) shows two periodic orbits in 
the configuration space. The motion at the upper end of 
the accessible energy interval is close to integrable. At 
a very high energy, motion in 0i direction is preferred, 
since the linear frequency wi of original mode 1 is higher 
than the frequency uj^ of mode 3. For decreasing energy, 
the KAM island around the center at -02 = increases in 
size while the one around ip2 = i^ decreases. The central 
periodic orbit around 02 = remains stable for ener- 
gies down to approximately i5 = 45, while the other one 
soon becomes unstable and its KAM island disappears. 
In Fig. 1(a) we see clearly the large KAM island belong- 
ing to the organization center 'ip2 = 0, with center at 
J2 = 4. In contrast to the idealized organization center 
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FIG. 2: (Color online) The classical reduced system for an 
energy E — 80. (a) Poincare plot as in Fig. 1(a). (b) Periodic 
orbits crossing the Poincare section in the centers of the KAM 
islands (black) for initial values {ipi,ip2, J2) = (0,0,1.3) and 
at the border (green) for initial values (0,7r,2.6). 

"02 = 0, the exact one is a periodic trajectory running 
in ipi showing small wiggles in -02 direction around the 
average value -02 — 0. However for our considerations 
it is simpler and completely satisfactory to replace this 
true organization center, the true periodic orbit, by an 
idealized organization center, for which we just take the 
straight line -02 = 0. The reader might remember the 
previous discussion of the onset of angle coupling and 
the values of the angles at which coupling sets in. In 
the spirit of this previous discussion we define the ideal- 
ized organization center as the subset of the configuration 
space defined by the angle restrictions exactly at the on- 
set of the corresponding coupling scheme. Also the ideal- 
ized semiclassical wave functions are given with respect 
to the corresponding idealized organization center. 
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Type (C) : If the effective frequencies of modes 1 and 
2 are close, but that of the third mode is not close, then 
we expect that modes 1 and 2 are locked but mode 3 
is independent. Then the motion is again quasiperiodic 
with two independent frequencies. In the reduced sys- 
tem, the organization center is a one dimensional curve 
which can be idealized by a line V'l = constant, where 
the constant usually is or tt according to the discussion 
in the beginning of this section. The periodic orbit it- 
self running in the '02 direction is almost impossible to 
find in Poincare maps with plane of intersection -01 = 0, 
since it violates the transversality of the map. However, 
when it is stable, then there is a bundle of invariant tori 
around it. In Poincare plots in the planes -01 = 0, these 
invariant tori appear as lines extending over all values of 
-02. In Fig. 1(a) they are the lines at the highest values 
of J2. In Fig. 3, we show some numerical results at en- 
ergy E = 40. Part (a) shows the Poincare map and parts 
(b) and (c) show trajectories in configuration space. In 
Fig. 3(a) the lines at small values of J2 belong to the 
tori around the organization center -01 = 0. Figure 3(b) 
shows a segment of a typical quasiperiodic orbit on one of 
these tori. While running monotonously in the negative 
-02 direction, it oscillates in -01 around the value 0. 
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FIG. 3: (Color online) The classical reduced system for an en- 
ergy E — 40. (a) Poincare plot as in Fig. 1(a). (b) Quasiperi- 
odic orbit with initial values {tpi,4'2,J2) = (0, tt, 2). (c) Two 
periodic orbits: one oscillating along the diagonal with start- 
ing point (0, 0, 14.4) (black) and the other rotating around 
the line 0i = 1/12 -I- ti" (green) with starting point (0, tt, 11.5). 

Type (D) : If the effective frequencies of modes 1 and 
3 are very close, then also the weak indirect tunneling 
processes between modes 1 and 3 can cause coupling. If 
the frequency of mode 2 is far from this common fre- 



quency, then mode 2 runs independently. The corre- 
sponding organization center in the reduced system is 
the line "01 = 0^2 + constant, where again this constant 
is usually or tt. In Poincare plots in planes "01 = 0, we 
see secondary islands. In Fig. 3(a), the two KAM island 
of moderate size with centers at ^2 — Q, J2 — 14 and 
-02 = TT and J2 = 11-5 respectively represent this type of 
motion. The two periodic orbits belonging to the centers 
of these two KAM islands are shown in Fig. 3(c). One 
oscillates along the diagonal and the other rotates around 
along the line 0i = ■02 + ti"- 

Type (E) : If all three effective frequencies are close, 
then there are two possibilities: 

(El): There is coupling between all three modes and 
the idealized organization center in the configuration 
space of the reduced system is a (fixed)point. The ac- 
tual trajectories oscillate around this coupling point and 
the relative angles V'fc do not rotate around the whole 
configuration torus. This behavior, which dominates at 
very small energy, is shown in Fig. 4 at energy E = 27. 
Only a limited range of "02 values around the point zero 
is energetically accessible. The same also holds for -01. 
Rotations around the configuration torus in either direc- 
tion or the diagonal become possible only for a higher 
energy. One of the organizing centers is represented in 
the Poincare plot by a stable fixed point which lies at 
the center of the large KAM island shown in Fig. 4(a), 
and which is shown in configuration space in Fig. 4(b) 
as the figure-of-eight orbit mainly oscillating in the an- 
tidiagonal direction. The other organizing center is an 
unstable periodic orbit belonging to the unstable fixed 
point near ^02 = 0, J2 = 12 in the Poincare plot. In 
the configuration space plot of Fig. 4(b), it is the orbit 
oscillating in the diagonal direction. At this energy, all 
trajectories in configuration space oscillate around the 
point (0,0), which acts as point organizing center. The 
two periodic orbits of Fig. 4(b) then act as guiding struc- 
tures for these fluctuations around the organizing center. 
Topologically speaking, all trajectories are contractible 
to a point on the configuration torus at very low energy. 
At the lower end of the accessible energy interval, the dy- 
namics starts as almost integrable and for this case the 
invariant manifolds of the unstable fixed point mentioned 
above lie close to a figure-of-eight shape separatrix in the 
Poincare section. For increasing energy the system moves 
further away from integrable and the separatrix breaks 
and turns into a homoclinic tangle which is the central 
structure of a chaos strip. In Fig. 4(a) for energy E = 21 
this chaotic layer still has moderate size. For higher en- 
ergy it grows rapidly and turns into the large chaotic sea 
seen in Fig. 3(a) at energy E = 40. 

(E2): The couplings break and reestablish intermit- 
tently, and the dynamics shows large-scale chaos. The 
appearance of chaos in the case of two independent reso- 
nant interactions becoming active is a demonstration of 
Chirikov's point of view of chaos being caused by reso- 
nance overlap [28]. In Poincare plots, we see large-scale 
chaos and eventually embedded in it remnants of islands 
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FIG. 4: (Color online) Classical reduced system for energy 
E — 27. (a) Poincare plot as in Fig. 1(a). (b) Trajectories 
through the stable fixed point (black, double loop) with initial 
values {tpi,tp2,J2) ~ (0,0,6.1), and the unstable fixed point 
(green, along the line tpi — ■02) with initial values (0,0, 12.1). 



and regular structures. The beginning of chaos for small 
energies can be seen in Fig. 4(a); chaos on a large scale 
is evident in Figs. 1(a) and 3(a). 



IV. CLASSIFICATION OF SEMICLASSICAL 
WAVE FUNCTIONS 

In this section, we show examples of wave functions be- 
longing to the various classes of motion described in the 
previous section. Our method of classification has been 
developed in [13-16] especially for Hamiltonians given 
quantum mechanically in raising and lowering operators 
and classically in action-angle variables. An analysis in a 
similar spirit of wave functions in the usual position space 
is rather common in molecular physics, two representa- 
tive examples are [29, 30]. In contrast to the procedure 
in action-angle space, the procedure in regular position 
space can be extended to scattering resonances, see [31]. 

Our strategy of classification is as follows: First we ex- 
pand the eigenstates |$) in the representation [■01, "02), 
according to Eq. (22). This representation of the eigen- 
states in the reduced configuration space is in complete 
analogy to the classical configuration space spanned by 
the angle variables '01 and ■02 and therefore allows a di- 
rect comparison between the classical and quantum sys- 
tem. We refer to these eigenfunctions $(0i,02) as the 
semiclassical wave functions in order to indicate this re- 
semblance. We then check whether the density of the 
semiclassical wave function in the reduced configuration 
space resembles the structure of one of the organization 
centers described in the previous section (types (A) - 
(E2)). I.e. we check, whether the density is distributed 
over the whole configuration space without clear nodal 
structures (type (A)), is concentrated along a few lines 
in the 0i direction (type (B)), in the 02 direction (type 
(C)) or in the diagonal direction (type (D)), is organized 
around the point center (0, 0) (type (El)) or shows ran- 
dom interferences between the pattern of different or- 



ganization centers leading to irregular structures (type 
(E2)). 

We call states, for which the density is located in a sin- 
gle crest along the organizing center, a transverse ground 
state to this organization center. In transversely excited 
states, the density is concentrated along various copies 
of the organization center, where these various copies are 
displaced relatively to each other and the wave function 
shows nodal structures between them. In addition, we 
look for the phase advance in directions in or parallel to 
the organization structure. The phase function must be 
continuous along curves which do not cross nodal lines. 
Recall that the phase function can have singularities only 
in zero points of the density. Accordingly, the curve along 
a crest of high density must be a curve of continuous 
phase. Then the phase advance of such a curve must be 
some integer multiple of 27r, say /i; • 27r, and this number 
fii serves as one quantum number of the state. These lon- 
gitudinal quantum numbers, together with the transverse 
quantum numbers given by the nodal structures, provide 
a complete set of quantum numbers characterizing the 
state relative to its organization center. We expect all 
states which can be related to an organization center to 
be close to a product of a plane wave in the longitudi- 
nal direction of this organization center and an oscillator 
function in transverse directions. 

In states belonging to classically chaotic motion, we 
do not see a simple and clear pattern in the density nor 
in the phase. Accordingly we are not able to give any 
assignment by quantum numbers to such states. 

In the following, the eigenstates ^k, k — 1, 2, . . .496 
are sorted by increasing energy starting with the label 1 
for the eigenstate with lowest energy. 



Point organization center (Type (El)) 

We start our analysis of the semiclassical wave func- 
tions at the lower end of the accessible energy interval. 
Since the Hamiltonian is dominated by quadratic anhar- 
monicities, the smallest energy is realized by distribut- 
ing the total excitation of 30 quanta (particles) evenly 
over the 3 basis modes (potential wells). In the classi- 
cal picture, this corresponds to the case where all three 
actions Ik are close to each other. Thus the three effec- 
tive frequencies (23) are very similar and frequency and 
phase locking is easily established by the resonant cou- 
pling terms in the Hamiltonian as explained in the begin- 
ning of section III. In the classical configuration space, 
this mechanism restricts the trajectories to a small region 
of configuration space (cf. Fig. 4) . This behavior is con- 
firmed in the quantum case. Here, the wave functions are 
organized around a point, as can be seen in Fig. 5, where 
the ground state and various excited states are plotted. 
State $1 is the ground state in this class. In this case 
the ground state of an organization center coincides with 
the energetic groundstate $i of the whole system, but 
we will assign also a groundstate for the other types of 
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guiding centers. The state $2 is the first transversal ex- 
citation in the antidiagonal direction, while the state $3 
represents the first transversal excitation in the diagonal 
direction. The state $4 represents the second transver- 
sal excitation in the antidiagonal direction, and the state 
$5 is the combination of one transverse excitation in the 
diagonal and one in the antidiagonal direction. State $9 
is the fourth excitation in the antidiagonal direction. A 
point center does not have any longitudinal directions. 
Accordingly, there are no phase advances in longitudinal 
directions to be counted for the assignment and any state 
of this class is characterized by the two transverse exci- 
tation numbers {ntd, l^ta) , one in the diagonal direction 
and one in the antidiagonal direction. Thus we show only 
the density plots without the phases in Fig. 5. 

In this scheme, the six states $1, $2, ^3, 'I'4, *&5 and $g 
have quantum numbers (0,0), (0,1), (1,0), (0,2), (1,1) 
and (0,4), respectively Note that the direction of exci- 



(a) 
1 


(b) 


(c) 
1 

•• 


(d) 


(e) 

A 


(f) 



\V^/n 



V^/Tt 



FIG. 5: Gray scale plot of the squared modulus of the (a) 
ground state <l?i('i/'i, V'2) and the excited states (b) $2, (c) 
"1>3, (d) $4, (e) $5 and (f) <l?g. White color corresponds to 
low density and black to the highest density. The range of 
the ipk is [-7r/2,37r/2]. 

tation corresponds to the direction of oscillation of the 
classical periodic orbits shown in Fig. 4(b). The classi- 
cal periods of these two orbits are Ta — 6.112 for the 
antidiagonal one and Td — 3.502 for the diagonal one. 



The quantum excitations in the corresponding direction 
increase the energy of the state by the classical frequency 
Lu = 27r/T, where T is the period of the orbit taken at an 
intermediate energy. 

The quantum-classical correspondence can be de- 
scribed in the following way: All three original modes 
are frequency locked and the phases fluctuate around the 
coupling point. The motion is similar to the one in a two 
dimensional anharmonic oscillator centered around the 
point (0, 0). This oscillator has its own normal modes and 
the states presented in Fig. 5 can be interpreted as some 
of the low lying excitation of this oscillator and described 
by the excitation numbers of these normal modes. How- 
ever, the reader should not confuse these modes of fluctu- 
ations around coupling points with the modes which are 
used to formulate the original Hamiltonian in Eqs. (3)- 
(5). Compare also with the discussion of the onset of 
coupling given in section III. 

The wave functions in this class are therefore close to 
two dimensional oscillator functions and can be described 
approximately by 

(26) 
where the functions Xn (x) are eigenfunctions of a one di- 
mensional oscillator with harmonic and anharmonic con- 
tributions. It is interesting to see what this means in the 
original coordinates ipk, Ik- Using the transformation 
(12), one obtains for the idealized eigenfunctions 

« e'""^' X,.. {Vi + ^3 - 2^2) Xm,„ {Vi - ^3) . (27) 

All three degrees of freedom arc entangled for this type of 
guiding center. The entanglement is the quantum analog 
of the phase locking in the classical picture. Altogether, 
we can assign 29 of the 496 eigenstates to this class of 
functions. 



Organization center i/>i ^ (Type (C)) 

The highest energies for a given number of particles 
are achieved by putting almost all excitation into one 
mode, with the other two modes having very low excita- 
tion. Classically, these two modes have similar effective 
frequencies, (see Eq. (23)), and therefore they are locked 
easily. In Fig. 6 we show as examples the densities and 
phases for the states $451 and $433. In part (a) we see 
the density concentrated along the line "01 = 0; thus, the 
transverse excitation number is /it — 0. Along this line, 
the phase function is almost like a plane wave. The total 
phase advance along one cycle around the organization 
center is 26 • 2tt. Accordingly, the longitudinal excitation 
number is fii = 26. Note that the phase function has 
singular points far away from the places of high density. 
In part (c) of the figure, we see the density concentrated 
along four lines in the 1/12 direction. The four density 
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crests are separated by 3 nodal lines, which can be seen 
very clearly as lines of discontinuities in the phase plot 
in part (d). Accordingly, the transverse excitation num- 
ber of state $433 is /it = 3. Along the density crests 
we count the total phase advance to obtain the longitu- 
dinal quantum number fii = 24. The energy distance 



The eigenfunctions in this class can therefore be writ- 
ten approximately as 




FIG. 6; Plot of the eigenfunctions $46i and $433 of the quan- 
tum system belonging to the V"! = guiding center. Plot (a) 
shows |<l?46ip, (b) shows arg($46i) mod27r, (c) shows |"l>433p, 
(d) arg($433) mod27r. In the phase plots, the degree of dark- 
ness from white to black indicates the phase advance from 
to 2ir. 



between two states which differ by one unit in /i; and 
that have the same transverse quantum number is given 
by the frequency of the classical organizing center, the 
periodic orbit (central fiber in the quasiperiodic motion 
shown in Fig. 3(b)), taken at an intermediate energy. 

The classical motion behind this class of states is as 
follows: Mode 3 runs with its own effective frequency 
independently of the other modes. The quantum number 
fii is its action due to the semiclassical assignment of the 
phase function ri{')pi,ip2) to the classical action integral. 



??(V'l,V'2) 



J ■ dip , 



(28) 



with J = (Ji, J2) and ip — [ipi, ip2). The phase ?7(V'i, V'2) 
is defined by $ = |$| exp{iri} and 7 is the classical guid- 
ing center. The path 7 is simply the line V'l = for this 
class and the the number of particles in mode 3 can be 
directly assigned to the quantum number /ii. Modes 1 
and 2 run locked with total excitation, i. e. number of 
particles, N — fii, and the quantum number fit charac- 
terizes the fluctuations of the coupled motion around the 
coupling point. 



$^„^.(V'i,^2)«e'^''e^^"^^XM.(^i). 



(29) 



Now we transform this expression back to the original co- 
ordinates, where we can interpret the actions Ik directly 
as the number of particles in the potential well k. Us- 
ing again transformation (12), we can write the idealized 
wave functions of this class as 



$,, 



Mi^l^^2,(P3) 



e^^'^^e^(^-^')^^XM*(^i-^2). (30) 



This type of wave function shows entanglement between 
modes 1 and 2, while mode 3 separates. The number of 
particles in mode 3 is given by /i/ while the transversal 
quantum number /it describes the transversal excitation 
of the organization center. A total of 51 eigenstates can 
be assigned to this class of functions. 



Organization center i/)2 = (Type (B)) 

The states of this class look very similar to those 
in the previous subsection, only with the roles of the 
modes 1 and 3 interchanged and hence with ipi and 
■!/'2 interchanged. However, there is no perfect symme- 
try between classes C and B because there is no perfect 
equality between the modes 1 and 3. Remember that 
LUi — —UJ3 y^ UJ3. This small perturbation of the sym- 
metry is responsible that the states of class B loose their 
characteristics under smaller transverse excitations as the 
ones for class C. Accordingly we can assign less states to 
class B, namely 42 only, than we have assigned to class 
C. 



Organization center ^j}^ = ^p2 (Type (D)) 

If almost all the action K is in mode 2, then modes 1 
and 3 have low actions and similar effective frequencies, 
whereas mode 2 has a quite different effective frequency. 
Even though the Hamiltonian does not contain a direct 
coupling between modes 1 and 3, sometimes the small 
indirect coupling is sufhcient to cause locking between 
modes 1 and 3. Fig. 7 shows the states $420 and $359 as 
two examples of semiclassical wave functions in this class. 
The organization center is the diagonal ipi = "02- State 
<i>42o has the transverse quantum number /it — relative 
to this center and state ^359 has fit — 1. The phase 
functions show that fii = 6 for state <&42o and fii — 8 
for state $420 • The energy distance between two states, 
which differ by one unit in /ij and that have the same 
transverse quantum number, is given by the frequency of 
the classical organizing center, namely the periodic orbit 
shown in Fig. 3(c) taken at an intermediate energy. 

The classical motion carrying these states is the fol- 
lowing: The coupled motion of modes 1 and 3 has the 
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number of particles /i/ while the rest of the total exci- 
tation N ^ III is in mode 2. The transverse quantum 
number fit again characterizes the fluctuations around 
the coupling point. For the idealized wave functions of 
the reduced system, we obtain 

1>M,,^.(^i, V'2) ~ e^"^" e'M,(^i+V2)/2 ^^^ (^^ _ ^^) (31) 

In the original coordinates, the wave function has the 
form 

In these coordinates, mode 2 separates from the other 
modes which are entangled. The number of particles in 
mode 2 is given by iV — /i;, while the rest of the parti- 
cles is in the entangled state of the other two modes, for 
which the quantum number fit is a measure of the fluc- 
tuations around the organization center. We can assign 
8 eigenstates to this class of functions. 
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FIG. 7: Plot of the eigenfunctions $420 and $359 of the quan- 
tum system belonging to the i/)! = ^2 guiding center. Plot 



(a) shows |$42o| , (b) arg(<J>42o) 
arg($359)mod27r. 



mod27r, (c) shows |$3 



(d) 



the two quantum numbers we count the phase advances 
around the two fundamental cycles of the toroidal con- 
figuration space. In part (b) of the figure we assign the 
quantum numbers fin = 2, fii2 = 5 and from part (d) we 
read off fin = 4 and fii2 = 1. 

These states are described by the classical motion in 
the following way: The original mode 1 has the number of 
particles fin and original mode 3 has fii2 particles. The 
rest of the excitation N — fin — fJ.12 is in mode 2. All 
three modes run independently with their own effective 
frequency. Thus phase functions of states of this class 
come close to a basis function (i.e. they resemble a plane 
wave), even though the wave function can be a strong 
mixture of several basis functions. The functional form 
of such states is therefore approximately given by 



$, 



i(V'l,V'2) 



^iN-ff gi(Mil'0l+Ml2'02) 



(33) 



or written in the original coordinates as 

$Mu.W2('^l, '/'2, <^3) - e^^"'"^ e'iN-^n-^i2W2 g.W2V.3 

(34) 
These idealized functions factorize and the three degrees 
of freedom are completely disentangled. There are 50 
eigenstates in this class. 
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FIG. 8: Plot of the eigenfunctions <I>4oi and $442 of the 
quantum system belonging to the T^ guiding center. Plot 
(a) shows |$40iP, (b) arg($40i) mod27r, (c) |<I>442p, (d) 
arg(<J>442) mod 2-k. 



Organization center T^ (Type (A)) 



Fig. 8 shows the wave functions of the states <i>4oi and 
$442, which do not show any coupling. These states be- 
long to normal mode motion in the original modes. This 
does not necessarily mean that they have a constant den- 
sity, but the density is without any clear structure and 
the phase function is close to a plane wave globally. As 



States based on chaotic motion (Type (E2)) 

Finally, we give two examples of wave functions where 
we could not make any assignment to one of the organiz- 
ing centers listed in the previous section. Fig. 9 shows 
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the densities and phases of states $100 and $i46- Neither 
in the density plots nor in the phase plots, can we dis- 
cover any clean pattern related to one of the organizing 
centers. The connection to the classical motion we inter- 
pret as follows: In classical chaos, any typical trajectory 
jumps around irregularly between the neighborhoods of 
various simple periodic orbits and therefore between var- 
ious types of motion. The corresponding quantum wave 
function should be random interference patterns of the 
structures belonging to the various organizing centers in- 
volved in the classical chaotic motion. Sometimes we can 
demix these interference patterns by forming appropri- 
ate linear combinations of several eigenfunctions of the 
Hamiltonian. 




FIG. 9: Plot of the eigenfunctions $100 and $146 of 
the quantum system belonging to class (E2). Plot (a) 



shows |$iooP, 



(b) arg($ioo)mod27r, (c) shows |$i46p, (d) 



arg(<I>i 



I mod 2iT. 



Concluding this section, we are able to characterize 
180 of the 496 eigenstates within the scheme of guiding 
centers given by the classical motion (excluding chaotic 
motion). Our aim is not a complete assignment of all 
states, but rather to give an easy visual criterion in order 
to select states with different types of e.g. entanglement 
and localization properties as described in this section 
for each class. For these states, one can use the classical 
picture in order to understand the quantum mechanical 
structure, which allows a very intuitive treatment of the 
states. The above graphical classification of the semi- 
classical wave functions is not strict and some functions 
allow ambiguous assignments. Such functions show char- 
acteristics of different classes and it is only a matter of 
degree in which class to put them. For example, the 
phase functions in Fig. 7 could be interpreted as continu- 
ous deformations of plane waves and therefore they could 



be assigned to type (A) as well. 



V. COMPARISON OF THE TIME DYNAMICS 

Finally, we wish to discuss the implications of our anal- 
ysis for the time evolution in the classical description 
of the system. The classical system can be interpreted 
as an array of three Bose-Einstein condensates where 
the condensate in each well is described by the Gross- 
Pitaevskii equation and where the condensates interact 
weekly through Josephson tunneling [3, 7]. 

In the previous section, we have used the classical sys- 
tem only to provide a tool for the classification of the 
quantum wave functions, and we have shown how close 
the quantum eigenfunctions resemble the classical guid- 
ing centers. In this section, we look in the other direc- 
tion. Starting from the classical system, i.e. the mean- 
field equations, we want to ask what information the 
structure of the quantum system can provide in order to 
solve the mean- field equations: the analysis of a system of 
coupled nonlinear differential equations is very involved, 
while in the quantum system we only have to diagonalize 
the Hamiltonian numerically and plot the eigenfunctions 
in configuration space. 

Since it is more convenient in this context to speak 
about complex occupation amplitudes, we introduce the 
new variables 



Cfc = Vh e 



(35) 



In these variables, the classical Hamiltonian (8) can be 
written as 



H 



E( 



LOklCkl 



- — (C1C2 



-Xk\Ck\ 



C2C1J 



^23. , 
— (C2C3 



C3C2J 



(36) 



with canonically conjugate variables (cfc,ic^) and corre- 
sponding equations of motion 



Cfc 



dH 



^Ck 



dH 

dck 



(37) 



This system of three ordinary differential equations for 
the complex coefficients Cfc is equivalent to the six equa- 
tions for the angles ipk and the actions Ik with k — 
1,2,3. The equations can also be derived from the Gross- 
Pitaevskii equation in coordinate space using an expan- 
sion of the condensate wave function in Wannier func- 
tions [32]. One can also use the eigenstates of the one- 
particle Hamiltonian, the so-called Wannier-Stark func- 
tions, resulting in the disappearance of the linear tun- 
neling terms in the Hamiltonian (3), while higher order 
coupling terms become important [11]. 

Now we choose the initial conditions Cfc (t = 0) by using 
the semiclassical correspondence (7) between the classical 
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actions Ik and the quantum numbers n^ of a number 
state 1^1, 7i2, 713): 



h 



1 



(38) 



In this way, we can construct initial conditions Cfe(O) — 
y/Tk, where the action Ik can be interpreted quantum 
mechanically via Eq. (38) as the number of particles in 
mode k. Furthermore, we can use this correspondence 
in order to construct initial conditions resembling the 
properties of the eigenstates of the system. Before we 
explain this in more detail we first discuss the case of the 
basis vectors. 



A. Basis vectors 

Here we investigate to which extend we can attribute 
the same characteristics to the quantum mechanical num- 
ber states |?T.i, 712, JT-a) and their classical analog defined 
by Eq. (38). Accordingly, we define the initial conditions 
for the time evolution of Eq. (37) as 



b{t = 0; ni, 712, ns) = (ci(0), 03(0), 03(0)) 



{^ni + l/2, V"2 + 1/2, V"3 + 1/2) . 



(39) 



In the following we will not explicitly write down the 
dependence of 6(i; tii, 712, 713) on the initial condition 
through the parameters (711,712,713) and simply use b{t). 
We fix the three initial phases to zero, which corresponds 
to zero imaginary part of the Cfc(O). With this initial con- 
ditions the time evolution can be calculated numerically, 
as shown in Fig. 10 for initial values (-\/2.5, \/5.5, \/23.5) 



C2 are locked, while C3 evolves independently. The differ- 
ence in the amplitudes between mode 3 and the other two 
prohibits a coupling. The amplitudes show a quite regu- 
lar oscillation in all three modes. This is motion of type 
(C) introduced in section III. Physically interpreted, the 
wells 1 and 2 couple through Josephson tunneling and 
the population between the two wells is exchanged peri- 
odically. In contrast, the number of particles of well 3 
stays approximately constant and much higher than the 
population of the other wells. This behavior reflects the 
well-known macroscopic self-trapping found in the dou- 
ble well potential [3] . Another type of this self-trapping 
effect in the type (C) dynamics can occur, when wells 1 
and 2 have approximately the same population iV/2 and 
well 3 is nearly empty. One can also observe the other 
types of dynamics in the vectors hit), except type (D), 
due to the very weak indirect coupling between modes 1 
and 3. The different time evolutions b{t) can be easily 
assigned to the different guiding centers by looking at the 
phases: 
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FIG. 10: (Color online) Time evolution of Eq. (37) for an ini- 
tial condition fo(0) = (\/2.5, \/5.5, \/23.5). Shown are squared 
modulus (top) and the phase of the first (solid, black), sec- 
ond (dashed, green) and third (dash-dotted, red) mode. In 
the phase plot the first and second phase almost coincide and 
lie above the third phase which has a bigger phase velocity. 
The time is measured with respect to T = 2-n /uj. 

using Eq. (37). In this example the phases of ci and 



FIG. 11: (Color online) Time evolution for an initial condition 
b{t = 0) = (V23.5, \/7.5, \/0.5). Shown are squared modulus 
(top) and the phase of the first (solid, black), second (dashed, 
green) and third (dash-dotted, red) mode. The time is mea- 
sured with respect to T = 2-k /uj. 

Type (A): All three phases behave independently and 
the amplitudes oscillate regularly. The individual con- 
densates in the different wells are completely decoupled 
and the population in each well stays approximately con- 
stant. 

Type (B): The dynamics shows the same behavior as 
for type (C), but with phase locking between mode 2 and 
3. 

Type (D): This type of motion is difficult to identify, 
because the indirect phase locking between modes 1 and 
3 is very weak. This leads to the effect that the phase 
velocities of these two phases are very close, but still dis- 
tinguishable. This is of course not a strict statement, and 
it depends on how long the time propagation is consid- 
ered. The problems with the classification of this type 
can also be seen in the quantum case in Fig. 7. In parts 
(b) resp. (d), the phase singularities are not sharp but 
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rather smooth, so these states could be assigned to type 
(A) as well. 

Type (El): In this case all three phases evolve with the 
same velocity and the amplitudes show similar regular 
oscillations as in types (B) and (C) for two locked phases. 

Type (E2): This class is characterized by intermitten- 
cies as illustrated in Fig. 11. The dynamics can be in- 
terpreted in such a way that the trajectories jump irreg- 
ularly between different coupling schemes. Accordingly, 
frequency locking between different pairs of modes is only 
established temporarily during the time evolution. 

With this scheme, we can classify the dynamics of all 
possible basis states 6(i), as shown in Fig. 12. The in- 
teresting point is that we can compare these results with 
the information that we extract from the semiclassical 
wave functions. For this we compare for a given basis 
state |ni, n2, n^) all eigenfunctions (22) to which the ba- 
sis state contributes significantly and assign a type (A)- 
(E2) to this basis state if possible. The result is shown in 
Fig. 13. The points with no symbol indicate states which 
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FIG. 12: 

actions Ik 



(Color online) Characterization of the classical 



\ck\ through direct numerical integration of 
Eq. (37). The action h is given hy h^ K-h-h- Plotted 
are time evolutions of type (A) (o, red), type (B) (V, green), 
type (C) (A, blue), type (El) (*, black) and type (E2) (□, 
cyan). 



cannot be assigned uniquely to a certain type. However, 
for the shown basis states one can see a close correspon- 
dence between the classical and the quantum system. 
Only at the fringes are there small deviations. There- 
fore the quantum mechanical analysis provides a grid of 
initial conditions for which we can predict the behav- 
ior of the solutions of the mean-field equations. Finally, 
we remark, that the classification of the basis states in 
Fig. 12 holds in principle also for an arbitrary choice of 



the initial phases in Eq. (39). Only at the fringes of the 
different zones does the behavior of the time dynamics 
depend crucially on the initial conditions and there it 
can deviate from this classification. 
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FIG. 13: (Color online) Characterization of the classical ac- 
tions Ik = \ck\ through the semiclassical wave functions. The 
action I2 is given by /2 = K — Ii — I3. Plotted are actions 
whose quantum analog belongs to type (A) (o, red), type (B) 
(V, green), type (C) (A, blue), type (D) (o, magenta), type 
(El) (*, black), and type (E2) (D, cyan). 



B. Eigenstates 

In the last section we discussed the close resemblance 
between the quantum and the classical picture by assign- 
ing the same characterization scheme with types (A)- 
(E2) to the basis functions and the solutions of the mean- 
field equations. In this section we want to investigate, 
whether also the eigenstates of the quantum system can 
be reinterpreted classically, i. e. if they can be used to 
identify the different types of dynamical behavior in the 
system of the three Bose-Einstein condensates weakly 
coupled by Josephson junctions. We construct the clas- 
sical analog of Eq. (21) by defining the set of vectors 

^(m, 712, na) = (ni + 1/2, n2 + 1/2, 713 + 1/2) , (40) 



which are related to the vectors b{t = 0) by Bk = 
6^.(0) (cf. Eq. (39)). However, note that the vectors 
5(711,712,713), like the vectors 6(i; 711,712, 713) of Eq. (39), 
do not form a basis of C^. In analogy to Eq. (20) one can 
write 



$(t = 0) 






-6(711,712,713) 



(41) 



ni+n2+n3=N 
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FIG. 14: (Color online) Time evolution of the mean-field 
equations for an initial condition corresponding to the first 
quantum eigenstate. Shown are squared modulus (top) and 

the phase of the first (— , black), second ( , green) and 

third (— ■ — , red) mode. The time is measured with respect 
to T = l-KJui. 



where the real- valued coefficients c„j^„2_„3 are taken from 
Eq. (21). In this naive approach, the vector $ can be in- 
terpreted as the quantum expectation value of the action 

/ (/fe = nfe -I- 1/2) in the quantum state |<i>). 



($|/aJ$) = 



E 



(Tifc + 1/2) , (42) 



where we have simply used the representation (20) of the 
eigenf unctions. The initial phases are chosen equal zero 
like in the case of the basis vectors. In order to use this 
vector $ as initial conditions for the mean- field equations, 
we must take the square root of each component, and to 
this end we define the new vector with components 
^fc = V^fc ■ These vectors are normalized as 



3 

E' 

k=\ 



E 



3 



Bi,^K. 



(43) 



ni+n2+"3=-'V 



k=\ 



where K — 31.5 ~ N ^ 3/2 is the classically conserved 
total action of Eq. (10). In the context of the Gross- 
Pitaevskii equation, the norm of the condensate wave 
function gives the number of particles in the condensate. 
We get the additional term of 3/2 for the number of par- 
ticles compared to the many-particle Hamiltonian (3), 
since we use the semiclassical correspondence of Eq. (7). 
For Bose-Einstein condensates with a number of parti- 
cles much larger than 30, one can ignore the term 1/2 in 
Eq. (7) and obtain the standard correspondence between 
the particle numbers. However, for TV = 30, semiclassical 
studies like the present work show that the identification 
(7) gives a much better agreement between classical and 
quantum mechanics. In order to obtain the normaliza- 
tion |0'p = 1, one simply has to set </> = (\)' \fN and 
replace the nonlinearities x\^ by xu = gjK . 



In Fig. 14, the time evolution for the initial condition 
$1 is shown. The time evolution shows approximately 
constant occupations |cfep (upper panel), and the three 
phases are locked. In the reduced system, this corre- 
sponds to a point in the neighborhood of a fixed point. 
For the parameter values chosen in this article, there does 
not exist an exact fixed point of the Hamiltonian flow of 
the reduced system, although this point serves as guid- 
ing center for the wave functions of type (El). In that 
sense the semiclassical wave functions behave very simi- 
larly in the neighborhood of a guiding center, while the 
solutions of the Gross-Pitaevskii equation are very sen- 
sitive to small deviations due to the nonlinearity of the 
time-evolution. 

Another example is shown in Fig. 15 for a type (A) 
motion. The phases of the modes evolve independently 
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FIG. 15: (Color online) Time evolution of Eq. (37) for the 
initial condition "1>444(0). Shown are squared modulus (top) 

and the phase of the first (— , black), second ( , green) and 

third (— ■ — , red) mode. The time is measured with respect 

to T = I-KJUJ. 



and the amplitudes show tiny oscillations, due to the fact 
that the time evolution does not coincide with the cor- 
responding idealized guiding center of type (A) . Because 
the system is dominated by the anharmonicities the effec- 
tive frequencies are almost linear in the actions according 
to Eq. (23). Therefore the slopes of the phase curves are 
proportional to the average values of the corresponding 
actions. 

To conclude, from the classical point of view the anal- 
ysis of the corresponding quantum system offers a direct 
visual method for the understanding of the structure and 
can be used to identify the dynamical behavior of the 
system of the three weakly coupled Bose-Einstein con- 
densates in the mean-field approximation simply by di- 
agonalizing the quantum Hamiltonian and plotting the 
eigenfunctions in the appropriate basis. 
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VI. CONCLUSION 

In our investigation of a Bosc-Einstcin condensate in 
a multi-well potential, we showed a close correspondence 
between the quantum mechanical description and a clas- 
sical version where the bosonic creation and annihilation 
operators of the many particle system are replaced by c- 
numbers. We truncated the many-particle Hamiltonian 
to a few relevant modes and obtained a system of three 
coupled anharmonic oscillators. Whether the truncation 
at a small number of modes is justified depends crucially 
on an appropriate choice of the expansion basis and on 
the external potential. In order to compare the quan- 
tum system with its classical counterpart, we introduced 
the concept of the semiclassical wave functions defined 
on the same toroidal configuration space as in the classi- 
cal system. This choice of the quantum mechanical rep- 
resentation allowed us to compare the quantum system 
directly with the classical system. In both cases, for the 
classical and the quantum system, we used the conserved 
particle number resp. total action to reduce the degrees 
of freedom to two. Classically, we can identify various 
geometric structures in phase space that are connected 
to different types of motion in the configuration space. 
These different types of motion belonging to the various 
guiding centers, are also found in the quantum mechan- 
ical wave functions. So we used these guiding centers 
firstly to sort a large number of wave function into these 
different classes, and secondly to assign uniquely geomet- 
ric quantum numbers to the wave functions within one 
class. In this geometric picture, the wave functions de- 



scribe the quantum excitations of the underlying classical 
dynamics. As an application, we can characterize the en- 
tanglement between the different modes and we can also 
determine the number of particles in each of the entan- 
gled modes using their associated quantum numbers. 

In the last part of this article we analyzed the signif- 
icance of the quantum mechanical classification of the 
wave functions for the classical dynamics. For this we 
studied classical trajectories which have initial conditions 
corresponding to quantum mechanical number states, or 
which correspond to the eigenstates directly. In both 
cases, we could obtain the characteristics of the semi- 
classical classification also from the classical trajectories, 
although the classical dynamics is much more sensitive 
to deviations from the idealized guiding centers. 

Concluding, wc showed that semiclassical wave func- 
tions provide an intuitive picture of the quantum me- 
chanical many-particle eigenfunctions, and allow a direct 
classification of the dynamics. 
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